
Job: betamc.std.prg (/home/gov/faculty/ppaolino/parr)
Sat Jun 16 17:37:20 2001, pid=14722, uid=1141(ppaolino)

cheng alg used if cheng=1, johnk otherwise, cheng=       1.0000000 
parm1 =        2.0000000       0.30000000 
parm2 =        2.0000000      -0.30000000 
p =        7.6844199 
q =        7.7719799 
true first difference (mean)=      0.29106368 
true first difference (variance)=   0.00014470353 

descriptive stats for standard beta parameters

       2.0079960      0.063189974 
      0.30431639      0.065902568 
       2.0079055      0.062316498 
     -0.29509473      0.065737692 
 
standard beta mse=     0.014109201 
 
alternative beta mse=     0.014109088 
 
ols mse=     0.014167523 
 
logit trans ols mse=     0.014188343 
 
het ols mse=     0.014167702 
 
logit trans het ols mse=     0.014186274 
 
log trans het ols mse=     0.016132657 
 
first differences for mean and variance for std beta (sim)
      0.29067198      0.011324166 
   2.1286626e-05     0.0017140506 
 
first differences for mean and variance for std beta (nosim)
      0.29071693      0.011326728 
   2.4059928e-05     0.0017109446 
 
first differences for mean and variance for alt beta (sim)
      0.28896822      0.011207565 
  -5.5607407e-05     0.0017097517 
 
first differences for mean and variance for alt beta (nosim)
      0.28903262      0.011195669 
  -5.9142760e-05     0.0017085747 
 
first differences for mean ols
      0.27562753     0.0097198683 
 
first differences for mean ols (logit trans)
      0.31097293      0.011960961 
 
first differences for mean and variance for het normal
      0.27576486     0.0097911615 
  -0.00013743577     0.0018917281 
 
first differences for mean and variance for logit trans.
      0.31070708      0.012045711 
   0.00017285556     0.0022411312 
 
first differences for mean ols (log trans)
      0.29584219      0.012367823 
 
efficiency of standard beta vs. alternative beta=       100.41080 
 
efficiency of beta vs. ols=       161.03010 
 
efficiency of beta vs. ols (log trans)=       205.03324 
 
efficiency of standard beta vs. het normal=       160.34332 
 
efficiency of standard beta vs. het normal (log trans)=       203.41433 
 
efficiency of standard beta vs. OLS (log trans)=       117.01084 
 
efficiency of standard beta vs. alternative beta=       100.32107 
 
efficiency of standard beta vs. het normal=       111.51330 
 
efficiency of standard beta vs. het normal (log trans)=       130.67351 
 
standard beta parms' overconfidence=
       101.42127 
       104.83233 
       100.05407 
       104.67781 
 
alternative beta parms' overconfidence=
       103.31591 
       99.201841 
       101.11796 
       100.36163 
 
ols parms' overconfidence=
       102.70991 
       90.980038 
 
ols (log. trans) parms' overconfidence=
       102.37168 
       106.87561 
 
ols (log trans) parms' overconfidence=
       99.703866 
       107.46846 
 
heteroskedastic normal parms' overconfidence=
       102.94461 
       91.724748 
       94.398581 
       93.682903 
 
heteroskedastic log. trans. parms' overconfidence=
       102.56790 
       107.66981 
       105.65899 
       118.39102 
 
standard beta's mean first difference overconfidence
       103.23041 
 
standard beta's variance first difference overconfidence
       103.71659 
 
alternative beta's mean first difference overconfidence
       99.346448 
 
alternative beta's variance first difference overconfidence
       99.531340 
 
check standard beta's percentage within 95% conf int

      0.94000000 
      0.94000000 
      0.94800000 
      0.94000000 
 
check alternative beta's percentage within 95% conf int

      0.95300000 
      0.94800000 
      0.94300000 
      0.95600000 
 
check ols percentage within 95% conf int

      0.95300000 
      0.96800000 
 
check logit trans/ols percentage within 95% conf int

      0.95800000 
      0.93100000 
 
check het ols percentage within 95% conf int

      0.95300000 
      0.96600000 
      0.95300000 
      0.96900000 
 
check logit trans/ols percentage within 95% conf int

      0.95700000 
      0.92700000 
      0.93500000 
      0.91300000 
 
check logit trans/ols percentage within 95% conf int

      0.95500000 
      0.92500000 
 
check percentage of sig parms (standard beta)
       1.0000000 
     0.051000000 
 
check percentage of sig parms (alternate beta)

     0.047000000 
       1.0000000 
       1.0000000 
     0.048000000 
 
check percentage of sig parms (basic ols)

       1.0000000 
       1.0000000 
 
check percentage of sig parms (logit trans ols)

     0.042000000 
       1.0000000 
 
check percentage of sig parms (het ols)

       1.0000000 
       1.0000000 
       1.0000000 
     0.031000000 
 
check percentage of sig parms (het ols - logit trans)

     0.042000000 
       1.0000000 
       1.0000000 
     0.088000000 
 
check percentage of sig parms (ols -- log tranform)

       1.0000000 
       1.0000000 
 
descriptive stats for y's for successful runs

      0.49691536 
      0.18219963 
     0.061573412 
      0.94198230 
descriptive stats for y's for failed runs
       0.0000000 
descriptive stats for x's for successful runs

    -0.020263478 
      0.99912922 
      -2.7494157 
       2.9375649 
descriptive stats for x's for failed runs
       0.0000000 
number of failed runs
       0.0000000 

code used to generate the above results:

library cml;
#include cml.ext;
CMLset;

fails=0;
crash=0;
ystatf=0;
xstatf=0;

/* load x matrix based upon number of observations */

x = ones(500,1)~rndn(500,1);
n = rows(x);

/* produce matrices for simulated first differences */

meanx=meanc(x); @ create vector of means for the indep vars @
stdx=stdc(x);   @ create vector of stdev for the indep vars @
minx=minc(x);
maxx=maxc(x);

/* start obtaining the first difference vectors (means) */

sw=zeros(1,cols(x)-1)|eye(cols(x)-1);
xl=meanx-stdx.*sw;           @ these are kXk-1 vectors @
xh=meanx+stdx.*sw;       
xmin=meanx.*(1-sw)+minx.*sw;  
xmax=meanx.*(1-sw)+maxx.*sw;

/* start creating Monte Carlo data */
k = cols(x);
parm1 = {2.0, 0.3};
parm2 = {2.0, -0.3};
cheng=1;
print "cheng alg used if cheng=1, johnk otherwise, cheng=";;cheng; 
p = exp(x*parm1);
q = exp(x*parm2);

/* determine true first differences */

/* get  alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alt=exp(xl'*parm1);                          @ nxk-1 matrix @
         blt=exp(xl'*parm2);
         aht=exp(xh'*parm1);
         bht=exp(xh'*parm2);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlt=alt./(alt+blt);                    @ nxk-1 matrix @
         meanht=aht./(aht+bht);
         varlt=(alt.*blt)./(((alt+blt)^2).*(alt+blt+1));
         varht=(aht.*bht)./(((aht+bht)^2).*(aht+bht+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmt=meanht-meanlt;
         fdvt=varht-varlt;

/* print data parameters for simulation */

print "parm1 = ";; parm1';
print "parm2 = ";; parm2';
print "p = ";; meanc(exp(x*parm1));
print "q = ";; meanc(exp(x*parm2));
print "true first difference (mean)=";; fdmt;
print "true first difference (variance)=";; fdvt;

/* start Monte Carlo trials */
trial=1;
do while trial<=1000; 

print "MC trial=";; trial;

   if trial==1;
      ret=0;
      if fails==0;
         bt=0; 
      endif;
   endif;

   if fails==0;
   success=0;

/* comment out random number generator if p,q<1 */
if cheng==1;
/* create J&K (Chen) alpha */
ralpha=p+q;

/* create J&K (Chen) beta */
p1=p.>1;
q1=q.>1;
rbeta=real((p1.*q1).*(sqrt((ralpha-2)./(2.*p.*q-ralpha))) +
      (1-p1.*q1).*(1/(minc((p'|q')))));

/* create J&K (Chen) gamma */
rgamma=p+1/rbeta;

/* step 1 in J&K (Chen) RNG */
u1=rndu(n,1);
u2=rndu(n,1);
v=rbeta.*ln(u1./(1-u1));
w=p.*exp(v);

/* step 2 in J&K (Chen) RNG */

done=0;
      do while done<n;
         d=(ralpha.*ln(ralpha./(q+w))+rgamma.*v-1.3862944).<ln((u1^2).*u2); 
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v=(1-d).*v + d.*rbeta.*ln(u1./(1-u1));
         w=(1-d).*w + d.*p.*exp(v);
         done=n-sumc(d);
      endo;

/* step 3 in J&K (Chen) RNG */
y=w./(q+w); 
endif;

/* discretize y 

y1=yc.>0;
y2=yc.>0.14;
y3=yc.>0.28;
y4=yc.>0.42;
y5=yc.>0.56;
y6=yc.>0.70;
y7=yc.>0.84;

y=(y1+y2+y3+y4+y5+y6+y7)/8; */

/* comment out random number generator if p or q>1 */
if cheng/=1;
/* step 1 in Johnk's RNG */

u1=rndu(n,1);
u2=rndu(n,1);

v1=u1^(1/p);
v2=u2^(1/q);

/* step 2 in Johnk's RNG */

w=v1+v2;

/* step 3 in Johnk's RNG */

done=0;
      do while done<n;
         d=w.>1;
         u1=(1-d).*u1 + d.*rndu(n,1);
         u2=(1-d).*u2 + d.*rndu(n,1);
         v1=(1-d).*v1 + d.*u1^(1/p);
         v2=(1-d).*v2 + d.*u2^(1/q);
         w=v1+v2;
         done=n-sumc(d);
      endo;

/* step 4 in Johnk's RNG */
y=v1./w;
endif;

/* cheat function to re-set Ys that equal 1 --- would choke llik function
   otherwise */
d=y.==1;
print "cheat=";;sumc(d);
y=(1-d).*y + d.*.99999;  


/* return the monte carlo data set */
    data = y~x[.,2:cols(x)];

   endif;

   if fails==0;

      do while success==0 and fails<=30;

/* maximum likelihood function for estimating parameters alpha and beta */
__title=" ** A Beta Model, Analytic Solution (Standard Estimation) **";

stval=1.0|ones(k-1,1)*0.5|1.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

       if fails>=1 and fails<=40;
          stval=(bt')/2;
          print "retrying run";
       endif;

print "mean y, max y, min y, std y";
meanc(y)~maxc(y)~minc(y)~stdc(y);

proc psi(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;

proc gradproc(theta,ds);
    local y,x,z,be,h,xb,p,q,per1,gr1,gr2,grad;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    per1=p+q;
    gr1= psi(p+q).*x.*p-psi(p).*x.*p+x.*p.*ln(y);
    gr2= psi(p+q).*x.*q-psi(q).*x.*q+x.*q.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);
    endp;

proc loglik(theta,ds);
    local y,x,z,be,h,xb,p,q,l1,l2,l3,l4,l5,llik;
		     
    y=ds[.,1];
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    be=theta[1:cols(x),1];                 @ parameters for alpha @
    h=trimr(theta,rows(be),0);            @ parameters for beta @
    p=exp(x*be);                 @ expression for alpha @
    q=exp(x*h);                  @ expression for beta @

    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradproc;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400;

   if fails.>=5 and fails.<=10;
     _cml_DirTol=.000001;
   endif;
   if fails.>10 and fails.<=15;
     _cml_DirTol=.0000001;
   endif;
   if fails.>15 and fails.<=20;
     _cml_DirTol=.00000001;
   endif;
   if fails.>20;
     _cml_DirTol=.00001;
   endif;

    {b,logl,g,vc,ret}=CMLprt(CML(data,0,&loglik,stval));

/* check that parameters are reasonable */
   bt=(b');

/* determine success or failure of the trial */
   cor=vc[1,1]/vc[1,1];
         if ret>=1 or cor/=1;
            fails=fails+1;
         elseif ret==0 and cor==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

   if fails==0;
      success=0; bta=0;
      do while success==0 and fails<=30;

CMLset;
/* Alternative beta estimation */
__title=" ** A Beta Model, Numerical Solution (Alternative Estimation) **";

stvala=0.0|ones(k-1,1)|2.0|ones(k-1,1)*0.4;

/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvala=(bta')/2;
         print "retrying run";
      endif;

proc gradpa(theta,ds);
    local y,one,x,z,be,h,xb,zh,zh1,e,per1,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    zh1=zh+1;
    e=xb./(1+xb);
    per1=(-x.*e+x.*(e^2)).*zh1;

    gr1=psi(e.*zh1+(1-e).*zh1-1).*(x.*e.*zh1-(e^2).*zh1.*x+per1)-
        psi(e.*zh1-e).*(x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x)-
        psi((1-e).*zh1-1+e).*(per1+x.*e-(e^2).*x)+
        (x.*e.*zh1-(e^2).*zh1.*x-x.*e+(e^2).*x).*ln(y)+
        (per1+x.*e-(e^2).*x).*ln(1-y);
    gr2=psi(e.*zh1+(1-e).*zh1-1).*(e.*z.*zh+(1-e).*z.*zh)-
        (psi(e.*zh1-e).*xb.*z.*zh)./(1+xb)-
        psi((1-e).*zh1-1+e).*(1-e).*z.*zh+xb.*z.*zh.*ln(y)./(1+xb)+
        (1-e).*z.*zh.*ln(1-y);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglika(theta,ds);
    local y,one,x,z,be,h,xb,zh,e,v,p,q,llik,ds1;
		     
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=exp(x*be);
    zh=exp(z*h);
    e=xb./(1+xb);
    v=e.*(1-e).*(1/(zh+1));
    p=((e^2).*(1-e))./v-e;
    q=(e.*(1-e)^2)./v-(1-e);
    if (sumc((p.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((q.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc((y.>0)))<n;
       llik=miss(0,0); print "run failed"; 
    elseif (sumc(((1-y).>0)))<n;
       llik=miss(0,0); print "run failed"; 
    else;
       llik=lngm(p+q)-lngm(p)-lngm(q)+(p-1).*ln(y)+
         (q-1).*ln(1-y);  @ log-likelihood @
    endif;
    retp(llik);endp;

_cml_GradProc=&gradpa;
_cml_GradCheck=1;
_cml_Options = { none }; 
_cml_MaxIters=400; 
    {ba,logl,g,vca,ret}=CMLprt(CML(data,0,&loglika,stvala));

/* check that parameters are reasonable */
   bta=(ba');

/* determine success or failure of the trial */
   cora=vca[1,1]/vca[1,1];
         if ret>=1 or cora/=1;
            fails=fails+1;
         elseif ret==0 and cora==1;
            success=1; fails=0;
         endif;
      endo;
   endif;

/* linear heteroskedastic estimation */

   if fails==0;
      success=0; bth=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvalh=inv(x'*x)*x'*y|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=10;
         stvalh=(bth')/2;
         print "retrying run";
      endif;  

proc gradph(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ds[.,1];	     
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikh(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ds[.,1];    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradph;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {blh,logl,g,vch,ret}=CMLprt(CML(data,0,&loglikh,stvalh));

/* check that Monte carlo parameters are reasonable */
bth=(blh');

/* determine success or failure of the trial */
   corh=vch[1,1]/vch[1,1];
       if ret>=1 or corh/=1;
          fails=fails+1;
          elseif ret==0 and corh==1;
          success=1; fails=0;
       endif;
    endo;
endif;

/* linear heteroskedastic estimation (with logit tranformation) */

   if fails==0;
      success=0; btl=0;
      do while success==0 and fails<=30;

CMLset;

__title=" ** A Linear-Normal Model, Numerical Solution **";

stvall=inv(x'*x)*x'*(ln(y./(1-y)))|(-0.5*ones(cols(x),1));
       
/* for failed runs, use starting values of crashed iteration */

      if fails>=1 and fails<=30;
         stvall=(btl')/2;
         print "retrying run";
      endif;  

proc gradpl(theta,ds);
    local y,one,x,z,be,h,xb,zh,gr1,gr2,grad;
    y=ln(ds[.,1]./(1-ds[.,1]));
    x=ones(rows(y),1)~ds[.,2:cols(ds)]; 
    z=x;
    be=theta[1:cols(x),1];
    h=trimr(theta,rows(be),0);      
    xb=x*be;
    zh=exp(z*h);
    gr1=(y-xb).*x./zh;
    gr2=-0.5*(x-(x.*(y-xb)^2)./zh);
    grad=gr1~gr2;
    retp(grad);endp;

proc loglikl(b,ds);                     
    local xb,x1,y,llik,b1,g1,s,s2,u,ds1,one;
    y=ln(ds[.,1]./(1-ds[.,1]));    /* dependent variable */
    one=ones(rows(ds),1);
    x1=one~ds[.,2:cols(ds)];
    b1=b[1:cols(x1),1];   /* partition off parms for E(y) */
    xb=x1*b1;
    s=trimr(b,rows(b1),0);      /* partition off parms for V(e) */
    g1=one~ds[.,2:cols(ds)];
    s2=exp(g1*s);
    u=y-xb;             
    llik=-0.5*(ln(s2)+u.*u./s2);  /* likelihood function */
    retp(llik);endp;

_cml_GradProc=&gradpl;
_cml_GradCheck=1;
_cml_Options = { none };
_cml_DirTol=0.000001; 
_cml_MaxIters=500;

    {bll,logl,g,vcl,ret}=CMLprt(CML(data,0,&loglikl,stvall));

/* check that Monte carlo parameters are reasonable */
btl=(bll');

/* determine success or failure of the trial */
   corl=vcl[1,1]/vcl[1,1];
      if ret>=1 or corl/=1;
         fails=fails+1;
      elseif ret==0 and corl==1;
         success=1; fails=0;
      endif;
   endo;
endif;

/* if trial has failed 30 times */

   if fails==30;
      print "run failed after 30 attempts";
      crash=crash+1;

      if ystatf==0;
         ystatf=(meanc(y)~stdc(y)~minc(y)~maxc(y));
      elseif ystatf/=0;
         ystatf=ystatf|(meanc(y)~stdc(y)~minc(y)~maxc(y));
      endif;

      if xstatf==0;
         xstatf=(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|maxc(x[.,2]));
      elseif xstatf/=0;
         xstatf=xstatf~(meanc(x[.,2])|stdc(x[.,2])|minc(x[.,2])|
                 maxc(x[.,2]));
      endif;
      fails=0;
      clear y;
   endif;

/* collect data for a successful trial */
   if success==1;

/* collect y stats for a successful trial */
      if trial==1;
         ystats=meanc(y)~stdc(y)~minc(y)~maxc(y);
         xstats=(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=b;
         sesum=sqrt(diag(vc));
         bsuma=ba;
         sesuma=sqrt(diag(vca));
         bsumh=blh;
         sesumh=sqrt(diag(vch));
         bsuml=bll;
         sesuml=sqrt(diag(vcl));
      endif;
      if trial>=2;
         ystats=ystats|(meanc(y)~stdc(y)~minc(y)~maxc(y));
         xstats=xstats~(meanc(x[.,2:cols(x)])|stdc(x[.,2:cols(x)])|
                minc(x[.,2:cols(x)])|maxc(x[.,2:cols(x)]));
         bsum=bsum~b;
         sesum=sesum~(sqrt(diag(vc)));
         bsuma=bsuma~ba;
         sesuma=sesuma~sqrt(diag(vca));
         bsumh=bsumh~blh;
         sesumh=sesumh~sqrt(diag(vch));
         bsuml=bsuml~bll;
         sesuml=sesuml~sqrt(diag(vcl));
      endif;

/* vectors of first difference values for the j^th independent var */
         bns=b[1:cols(x),.];                @ nxk matrix @
         hns=b[cols(x)+1:rows(b),.];   @ nxk matrix @

/* calculate model mse */

         ahat=exp(x*bns);
	 bhat=exp(x*hns);
         yhatns=ahat./(ahat+bhat);
	 bnsmse=sumc((y-yhatns)^2)/n;
	 
/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         alns=exp(xl'*bns);                          @ nxk-1 matrix @
         blns=exp(xl'*hns);
         ahns=exp(xh'*bns);
         bhns=exp(xh'*hns);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlns=alns./(alns+blns);                    @ nxk-1 matrix @
         meanhns=ahns./(ahns+bhns);
         varlns=(alns.*blns)./(((alns+blns)^2).*(alns+blns+1));
         varhns=(ahns.*bhns)./(((ahns+bhns)^2).*(ahns+bhns+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmns=meanhns-meanlns;
         fdvns=varhns-varlns;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmns=fdmns;  @ k-1X1 matrix for mean @
            totvns=fdvns;  @ k-1X1 matrix for variance @
	    totbmse=bnsmse;
         endif;

         if trial>=2 and ret==0; 
            totmns=totmns~fdmns;  @ k-1X(trial) matrix for mean @
            totvns=totvns~fdvns;  @ k-1X(trial) matrix for variance @
	    totbmse=totbmse|bnsmse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasim=(rndmn(b,vc,1000));  @ nxk matrix, where n=# of draws @
         be=betasim[.,1:cols(x)];                @ nxk matrix @
         h=betasim[.,cols(x)+1:cols(betasim)];   @ nxk matrix @

/* get predicted alpha and beta values with the j^th independent var 
   at +/- 1 sd and min/max values, while other x's held at their means */
         al=exp(be*xl);                          @ nxk-1 matrix @
         bl=exp(h*xl);
         ah=exp(be*xh);
         bh=exp(h*xh);

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlow=al./(al+bl);                    @ nxk-1 matrix @
         meanhigh=ah./(ah+bh);
         varlow=(al.*bl)./(((al+bl)^2).*(al+bl+1));
         varhigh=(ah.*bh)./(((ah+bh)^2).*(ah+bh+1));

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdm=meanhigh-meanlow;
         fdv=varhigh-varlow;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totm=meanc(fdm);  @ k-1X1 matrix for mean @
            totv=meanc(fdv);  @ k-1X1 matrix for variance @
            totms=stdc(fdm);  @ k-1X1 matrix for mean @
            totvs=stdc(fdv);  @ k-1X1 matrix for variance @
         endif;

         if trial>=2 and ret==0; 
            totm=totm~meanc(fdm);  @ k-1X(trial) matrix for mean @
            totv=totv~meanc(fdv);  @ k-1X(trial) matrix for variance @
            totms=totms~stdc(fdm);  @ k-1X(trial) matrix for mean @
            totvs=totvs~stdc(fdv);  @ k-1X(trial) matrix for variance @
         endif;

/* get relevant stats for alternative method of beta estimation */
/* vectors of first difference values for the j^th independent var */
         bnsa=ba[1:cols(x),.];                @ nxk matrix @
         hnsa=ba[cols(x)+1:rows(ba),.];   @ nxk matrix @

/* calculate mse for alternative beta estimation */
  
         yhatnsa=exp(x*bnsa)./(1+exp(x*bnsa));
 	 bnsamse=sumc((y-yhatnsa)^2)/n;
	 
/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanlnsa=(exp(xl'*bnsa)./(1+exp(xl'*bnsa)))';   @ nxk-1 matrix @
         meanhnsa=(exp(xh'*bnsa)./(1+exp(xh'*bnsa)))';
         dlsa=(1/(exp(xl'*hnsa)+1))';
         dhsa=(1/(exp(xh'*hnsa)+1))';
         varlnsa=meanlnsa.*(1-meanlnsa).*dlsa;
         varhnsa=meanhnsa.*(1-meanhnsa).*dhsa;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdmnsa=meanhnsa-meanlnsa;
         fdvnsa=varhnsa-varlnsa;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totmnsa=fdmnsa;  @ k-1X1 matrix for mean @
            totvnsa=fdvnsa;  @ k-1X1 matrix for variance @
	    totbamse=bnsamse;
         endif;

         if trial>=2 and ret==0; 
            totmnsa=totmnsa~fdmnsa;  @ k-1X(trial) matrix for mean @
            totvnsa=totvnsa~fdvnsa;  @ k-1X(trial) matrix for var @
	    totbamse=totbamse|bnsamse;
         endif;

/* vectors of first difference values for the j^th independent var */
         betasima=(rndmn(ba,vca,1000));  @ nxk matrix, where n=# of draws @
         bea=betasima[.,1:cols(x)];                @ nxk matrix @
         ha=betasima[.,cols(x)+1:cols(betasima)];   @ nxk matrix @

/* get predicted values of mean and variance with the j^th independent var
   at +/- sd and min/max values, while other x's held at their means */
         meanla=(exp(bea*xl)./(1+exp(bea*xl)));   @ nxk-1 matrix @
         meanha=(exp(bea*xh)./(1+exp(bea*xh)));
         dl=(1/(exp(ha*xl)+1));
         dh=(1/(exp(ha*xh)+1));
         varla=meanla.*(1-meanla).*dl;
         varha=meanha.*(1-meanha).*dh;

/* produce first differences for the predicted mean and variance with 
   the j^th indep var at +/- sd, while other x's held at their means */
         fdma=meanha-meanla;
         fdva=varha-varla;

/* produce vector for average first differences */

         if trial==1 and ret==0; 
            totma=meanc(fdma);  @ k-1X1 matrix for mean @
            totva=meanc(fdva);  @ k-1X1 matrix for variance @
            totmsa=stdc(fdma);
            totvsa=stdc(fdva);
         endif;

         if trial>=2 and ret==0; 
            totma=totma~meanc(fdma);  @ k-1X(trial) matrix for mean @
            totva=totva~meanc(fdva);  @ k-1X(trial) matrix for variance @
            totmsa=totmsa~stdc(fdma);
            totvsa=totvsa~stdc(fdva);
         endif;

/* set parameter vectors (linear heteroskeadastic) */
         bhh=blh[1:cols(x),1];
         hh=trimr(blh,rows(bhh),0);

/* calculate mse */

         yhatlh=x*bhh;
	 lhmse=sumc((y-yhatlh)^2)/n;
	 
/* get predicted means */
         elh=xl'*bhh;
         ehh=xh'*bhh;

/* get predicted variances */
         vlh=exp(xl'*hh);
         vhh=exp(xh'*hh);

/* produce first differences */
         if trial==1 and ret==0;
         totfdmh=ehh-elh;
         totfdvh=vhh-vlh;
	 totlhmse=lhmse;
         elseif trial>=2 and ret==0;
         totfdmh=totfdmh~(ehh-elh);
         totfdvh=totfdvh~(vhh-vlh);
	 totlhmse=totlhmse|lhmse;
         endif;

/* set parameter vectors (linear logit tranformation) */
         blo=bll[1:cols(x),1];
         hl=trimr(bll,rows(blo),0);
	 
/* calculate mse for linear het log trans */
  
         yhatlhlt=exp(x*blo)./(1+exp(x*blo));
	 lhltmse=sumc((y-yhatlhlt)^2)/n;

/* get predicted means */
         eol=exp(xl'*blo)./(1+exp(xl'*blo));
         eoh=exp(xh'*blo)./(1+exp(xh'*blo));

/* get predicted variances */
         vol=((exp(xl'*blo)./(1+exp(xl'*blo))^2)^2).*exp(xl'*hl);
         voh=((exp(xh'*blo)./(1+exp(xh'*blo))^2)^2).*exp(xh'*hl);

/* produce first differences */
         fdml=eoh-eol;
         fdvl=voh-vol;

/* produce first differences */
         if trial==1 and ret==0;
         totfdml=fdml;
         totfdvl=fdvl;
	 totthmse=lhltmse;
         elseif trial>=2 and ret==0;
         totfdml=totfdml~fdml;
         totfdvl=totfdvl~fdvl;
	 totthmse=totthmse|lhltmse;
         endif;

/* run ols */

{ vnam,m,bols,stb,vcols,stderr,sigma,cx,rsq,resid,dbw } = ols(0,y,x);

/* calculate mse */

      yhat=x*bols;
      olsmse=sumc((y-yhat)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsum=bols;
         selsum=stderr;
      endif;
      if trial>=2;
         blsum=blsum~bols;
         selsum=selsum~stderr;
      endif;

/* obtain the first differences */

olsl=xl'*bols;
olsh=xh'*bols;
fdols=olsh-olsl;

      if trial==1;
         fdolsum=fdols;
	 totomse=olsmse;
      endif;
      if trial>=2;
         fdolsum=fdolsum~fdols;
	 totomse=totomse|olsmse;
      endif;

/* run ols (with logit trans) */

yl=ln(y./(1-y));

{ vnam,m,bolsl,stb,vcols,stderrl,sigma,cx,rsq,resid,dbw } = ols(0,yl,x);

/* calculate mse */

      yhatlt=exp(x*bolsl)./(1+exp(x*bolsl));
      ltmse=sumc((y-yhatlt)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumo=bolsl;
         selsumo=stderrl;
      endif;
      if trial>=2;
         blsumo=blsumo~bolsl;
         selsumo=selsumo~stderrl;
      endif;

/* obtain the first differences */

olslo=exp(xl'*bolsl)./(1+exp(xl'*bolsl));
olsho=exp(xh'*bolsl)./(1+exp(xh'*bolsl));
fdolsl=olsho-olslo;

      if trial==1;
         fdolsuml=fdolsl;
	 tototmse=ltmse;
      endif;
      if trial>=2;
         fdolsuml=fdolsuml~fdolsl;
	 tototmse=tototmse|ltmse;
      endif;

/* run ols (with log trans) */

ly=ln(y+.01);

{ vnam,m,bolsll,stb,vcols,stderrll,sigma,cx,rsq,resid,dbw } = ols(0,ly,x);

/* calculate mse */

      yhatll=exp(x*bolsll)-0.01;
      ltmsel=sumc((y-yhatll)^2)/n;

/* obtain the overconfidence scores for the OLS parameters */

      if trial==1;
         blsumol=bolsll;
         selsumol=stderrll;
      endif;
      if trial>=2;
         blsumol=blsumol~bolsll;
         selsumol=selsumol~stderrll;
      endif;

/* obtain the first differences */

olslol=exp(xl'*bolsll)-0.01;
olshol=exp(xh'*bolsll)-0.01;
fdolsll=olshol-olslol;

      if trial==1;
         fdsumll=fdolsll;
	 totmsell=ltmsel;
      endif;
      if trial>=2;
         fdsumll=fdsumll~fdolsll;
	 totmsell=totmsell|ltmsel;
      endif;

/* reset fails and add a trial for a successful trial */
      fails=0;
      trial=trial+1; 
      clear y;

   endif;

endo;

/* check that parameters are accurate */
print "descriptive stats for standard beta parameters";
meanc(bsum')~stdc(bsum'); 
print " ";
/* check the mse for the std beta */
  
print "standard beta mse=";;meanc(totbmse);
print " ";
print "alternative beta mse=";;meanc(totbamse);
print " ";
print "ols mse=";;meanc(totomse);
print " ";
print "logit trans ols mse=";;meanc(tototmse);
print " ";
print "het ols mse=";;meanc(totlhmse);
print " ";
print "logit trans het ols mse=";;meanc(totthmse);
print " ";
print "log trans het ols mse=";;meanc(totmsell);
print " ";

/* check mean and std. dev of first differences */

print"first differences for mean and variance for std beta (sim)";
meanc(totm')~stdc(totm');
meanc(totv')~stdc(totv');
print " ";

print "first differences for mean and variance for std beta (nosim)";
meanc(totmns')~stdc(totmns');
meanc(totvns')~stdc(totvns');
print " ";

print "first differences for mean and variance for alt beta (sim)";
meanc(totma')~stdc(totma');
meanc(totva')~stdc(totva');
print " ";

print "first differences for mean and variance for alt beta (nosim)";
meanc(totmnsa')~stdc(totmnsa');
meanc(totvnsa')~stdc(totvnsa');
print " ";

/* check for mean and std. dev. of first difference for OLS */
print"first differences for mean ols";
meanc(fdolsum')~stdc(fdolsum');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (logit trans)";
meanc(fdolsuml')~stdc(fdolsuml');
print " ";

print "first differences for mean and variance for het normal";
meanc(totfdmh')~stdc(totfdmh');
meanc(totfdvh')~stdc(totfdvh');
print " ";

print "first differences for mean and variance for logit trans.";
meanc(totfdml')~stdc(totfdml');
meanc(totfdvl')~stdc(totfdvl');
print " ";

/* check for mean and std. dev. of first difference for OLS (log trans) */
print"first differences for mean ols (log trans)";
meanc(fdsumll')~stdc(fdsumll');
print " ";

/* efficiency of expected values */
/* compare efficiency of standard beta with alternative beta */

effnba=((totmnsa-fdmt)^2);
effdba=((totmns-fdmt)^2);
effalt=100*(sqrt(sumc(effnba')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. alternative beta=";;effalt;
print " ";

/* compare efficiency of standard beta with ols */

effno=((fdolsum-fdmt)^2);
effdb=((totmns-fdmt)^2);
effols=100*(sqrt(sumc(effno')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols=";;effols;
print " ";

/* compare efficiency of standard beta with ols */

effnl=((fdolsuml-fdmt)^2);
effdb=((totmns-fdmt)^2);
effolsl=100*(sqrt(sumc(effnl')))./(sqrt(sumc(effdb')));
print "efficiency of beta vs. ols (log trans)=";;effolsl;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbh=((totfdmh-fdmt)^2);
effdba=((totmns-fdmt)^2);
effhet=100*(sqrt(sumc(effnbh')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. het normal=";;effhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

effnbll=((totfdml-fdmt)^2);
effdbal=((totmns-fdmt)^2);
efflog=100*(sqrt(sumc(effnbll')))./(sqrt(sumc(effdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;efflog;
print " ";

/* compare efficiency of standard beta with log trans normal */

effnbl=((fdsumll-fdmt)^2);
effdba=((totmns-fdmt)^2);
efflogl=100*(sqrt(sumc(effnbl')))./(sqrt(sumc(effdba')));
print "efficiency of standard beta vs. OLS (log trans)=";;efflogl;
print " ";

/* efficiency for variance term */
/* compare efficiency of standard beta with alternative beta */

veffnba=((totvnsa-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffalt=100*(sqrt(sumc(veffnba')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. alternative beta=";;veffalt;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbh=((totfdvh-fdvt)^2);
veffdba=((totvns-fdvt)^2);
veffhet=100*(sqrt(sumc(veffnbh')))./(sqrt(sumc(veffdba')));
print "efficiency of standard beta vs. het normal=";;veffhet;
print " ";

/* compare efficiency of standard beta with heteroskedasic normal */

veffnbll=((totfdvl-fdvt)^2);
veffdbal=((totvns-fdvt)^2);
vefflog=100*(sqrt(sumc(veffnbll')))./(sqrt(sumc(veffdbal')));
print "efficiency of standard beta vs. het normal (log trans)=";;vefflog;
print " ";

/* check standard beta parameters' overconfidence */

npoca=(bsum-meanc(bsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesum^2;
dpoc=sqrt(sumc(dpoca'));
ocp=100*(npoc./dpoc);
print "standard beta parms' overconfidence=";;ocp;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check alternative beta parameters' overconfidence */

npoca=(bsuma-meanc(bsuma'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuma^2;
dpoc=sqrt(sumc(dpoca'));
ocpa=100*(npoc./dpoc);
print "alternative beta parms' overconfidence=";;ocpa;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols parameters' overconfidence */

npoca=(blsum-meanc(blsum'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsum^2;
dpoc=sqrt(sumc(dpoca'));
oco=100*(npoc./dpoc);
print "ols parms' overconfidence=";;oco;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log. trans) parameters' overconfidence */

npoca=(blsumo-meanc(blsumo'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumo^2;
dpoc=sqrt(sumc(dpoca'));
ocol=100*(npoc./dpoc);
print "ols (log. trans) parms' overconfidence=";;ocol;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check ols (log trans) parameters' overconfidence */

npoca=(blsumol-meanc(blsumol'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=selsumol^2;
dpoc=sqrt(sumc(dpoca'));
ocoll=100*(npoc./dpoc);
print "ols (log trans) parms' overconfidence=";;ocoll;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic normal parameters' overconfidence */

npoca=(bsumh-meanc(bsumh'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesumh^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic normal parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check heteroskedastic logit transform parameters' overconfidence */

npoca=(bsuml-meanc(bsuml'))^2;
npoc=sqrt(sumc(npoca'));
dpoca=sesuml^2;
dpoc=sqrt(sumc(dpoca'));
ocph=100*(npoc./dpoc);
print "heteroskedastic log. trans. parms' overconfidence=";;ocph;
clear npoca,npoc,dpoca,dpoc;
print " ";

/* check overconfidence for the first differences (mean) */

print "standard beta's mean first difference overconfidence";
nmoca=(totm-meanc(totm'))^2;
nmoc=sqrt(sumc(nmoca'));
dmoca=totms^2;
dmoc=sqrt(sumc(dmoca'));
ocm=100*(nmoc./dmoc);
ocm;
print " ";

/* check overconfidence for the first differences (variance) */

print "standard beta's variance first difference overconfidence";
nvoca=(totv-meanc(totv'))^2;
nvoc=sqrt(sumc(nvoca'));
dvoca=totvs^2;
dvoc=sqrt(sumc(dvoca'));
ocv=100*(nvoc./dvoc);
ocv;
print " ";

/* check overconfidence for the first differences (mean) */

print "alternative beta's mean first difference overconfidence";
nmocaa=(totma-meanc(totma'))^2;
nmoca=sqrt(sumc(nmocaa'));
dmocaa=totmsa^2;
dmoca=sqrt(sumc(dmocaa'));
ocma=100*(nmoca./dmoca);
ocma;
print " ";

/* check overconfidence for the first differences (variance) */

print "alternative beta's variance first difference overconfidence";
nvocaa=(totva-meanc(totva'))^2;
nvoca=sqrt(sumc(nvocaa'));
dvocaa=totvsa^2;
dvoca=sqrt(sumc(dvocaa'));
ocva=100*(nvoca./dvoca);
ocva;
print " ";

/* check percentage of parameters within 95% conf int. */

print "check standard beta's percentage within 95% conf int";
cib=abs(bsum-meanc(bsum')).<(1.96*sesum);
sumc(cib')/(trial-1);
print " ";

print "check alternative beta's percentage within 95% conf int";
ciba=abs(bsuma-meanc(bsuma')).<(1.96*sesuma);
sumc(ciba')/(trial-1);
print " ";

print "check ols percentage within 95% conf int";
cio=abs(blsum-meanc(blsum')).<(1.96*selsum);
sumc(cio')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cit=abs(blsumo-meanc(blsumo')).<(1.96*selsumo);
sumc(cit')/(trial-1);
print " ";

print "check het ols percentage within 95% conf int";
cih=abs(bsumh-meanc(bsumh')).<(1.96*sesumh);
sumc(cih')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cil=abs(bsuml-meanc(bsuml')).<(1.96*sesuml);
sumc(cil')/(trial-1);
print " ";

print "check logit trans/ols percentage within 95% conf int";
cill=abs(blsumol-meanc(blsumol')).<(1.96*selsumol);
sumc(cill')/(trial-1);
print " ";

/* get information about percentage of significant parameters */
/* mean values */
print "check percentage of sig parms (standard beta)";
tstatbm=(abs(totm./totms)).>1.96;
tstatbv=(abs(totv./totvs)).>1.96;
sumc(tstatbm')/(trial-1);
sumc(tstatbv')/(trial-1);
print " ";

print "check percentage of sig parms (alternate beta)";
tstatba=(abs(bsuma./sesuma)).>1.96;
sumc(tstatba')/(trial-1);
print " ";

print "check percentage of sig parms (basic ols)";
tstatbo=(abs(blsum./selsum)).>1.96;
sumc(tstatbo')/(trial-1);
print " ";

print "check percentage of sig parms (logit trans ols)";
tstatbl=(abs(blsumo./selsumo)).>1.96;
sumc(tstatbl')/(trial-1);
print " ";

print "check percentage of sig parms (het ols)";
tstatboh=(abs(bsumh./sesumh)).>1.96;
sumc(tstatboh')/(trial-1);
print " ";

print "check percentage of sig parms (het ols - logit trans)";
tstatblh=(abs(bsuml./sesuml)).>1.96;
sumc(tstatblh')/(trial-1);
print " ";

print "check percentage of sig parms (ols -- log tranform)";
tstatbm=(abs(blsumol./selsumol)).>1.96;
sumc(tstatbm')/(trial-1);
print " ";


/* obtain information about y's */
print "descriptive stats for y's for successful runs";
meanc(ystats);
print "descriptive stats for y's for failed runs";
meanc(ystatf);

/* obtain information about x's */
print "descriptive stats for x's for successful runs";
meanc(xstats');
print "descriptive stats for x's for failed runs";
meanc(xstatf');

print "number of failed runs";
crash;
output off;